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The ensemble Kalman filter (EnKF) ( Evensenj 2009a I has proven effective in quantifying uncertainty in a 
number of challenging dynamic, state estimation, or data assimilation, problems such as weather forecasting 
and ocean modeling. In these problems a high-dimensional state parameter is successively updated based 
on recurring physical observations, with the aid of a computationally demanding forward model that prop- 
agates the state from one time step to the next. More recently, the EnKF has proven effective in history 
matching in the petroleum engineering community (Evensen 2009b Oliver and Chen 2010). Such applica- 
tions typically involve estimating large numbers of parameters, describing an oil reservoir, using data from 
production history that accumulate over time. Such history matching problems are especially challenging 
examples of computer model calibration since they involve a large number of model parameters as well as a 
computationally demanding forward model. More generally, computer model calibration combines physical 
observations with a computational model ~ a computer model - to estimate unknown parameters in the 
computer model. This paper explores how the EnKF can be used in computer model calibration problems, 
comparing it to other more common approaches, considering applications in climate and cosmology. 

Keywords: computer experiments; model validation; data assimilation; uncertainty quantification; Gaussian 
process; parameter estimation; Bayesian statistics 



1 Introduction 



The ensemble Kalman filter (EnKF) has proven effective in quantifying uncertainty in a number 
of challenging dynamic, state estimation, or data assimilation, problems. Applications include 
weather forecasting ( Houtekamer et al. , 2005 ) , ocean modeling ( Evensen , 2003 ) , storm tracking 



(Aksoy et al. , 2009), hydrology (Moradkhani et al. , 2005) and wildfire modeling (Mandel et al 



2004), just to name a few. In these data assimilation problems, a high-dimensional state parameter 



is successively updated based on recurring physical observations, with the aid of a computationally 
demanding forward model that propagates the state from one time step to the next. The EnKF 
iteratively updates an ensemble of state vectors, using a scheme motivated by the standard Kalman 



filter (Meinhold and Singpurwalla , 1983 West and Harrison, 1997), producing an updated ensemble 



of states that is affected by both the forward model and the physical observations. More recently, the 



EnKF has proven effective in history matching in the petroleum engineering community (Evensen 



2009b; Oliver and Chen, 2010). Such applications typically involve estimating large numbers of 



parameters, describing an oil reservoir, using data from production history that accumulate over 
time. Such history matching problems are especially challenging examples of computer model 
calibration since they involve a large number of model parameters as well as a computationally 
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demanding forward model. Unlike standard data assimilation problems, here focus is on estimation 
of a static model parameter vector, rather than an evolving state vector. 

This paper explores how the EnKF can be used in computer model calibration problems that 
typically have a collection of model parameters to be constrained using physical observations. 
We first use a simple 1-d inverse problem to describe standard Bayesian approaches to produce 
a posterior distribution for the unknown model parameter vector, as well as the resulting model 
prediction. We then go on to describe how the EnKF can be used to address this basic problem, with 
examples taken from the literature in climate and cosmology. We end with conclusions summarizing 
the strengths and weaknesses of using the EnKF for computer model calibration. 

1.1 A simple inverse problem 



A simple inverse problem in which rj{-), given 
by the black line, denotes the forward model, 
mapping the unknown parameter 9 into an ob- 
servable rj{6). 

The physical observation j/ is a noisy version of 

where e ~ N{0,ay = .1^). The horizontal gray 
line and band denote the measurement {y = .8) 
and its uncertainty y±2ay. 

The model parameter 6 is given a A'^(0, 1) prior. 
The resulting posterior for is given by the 
shaded density at the bottom of the figure. 



Figure 1: A simple 1-dimensional inverse problem and resulting posterior density. 

In order to describe the basic approaches to inverse problems, we first describe a simple, 1- 
dimensional inverse problem shown in Figure [!} We take r]{-) to denote the forward model. It 
requires a single model parameter 9, producing a univariate output rj(6) that is comparable to a 
physical measurement y. Here we take the sampling model for y to be normally distributed about 
the forward model's output when the true value of the model parameter is input 

L{y\rjie))^exp{-^a-\y-viO)f} 

where the observation error is assumed known to be ctj^ = .1. 

After specifying a standard normal prior for the model parameter, with tt{6) denoting the prior 
density, the posterior density is given by 

Trie\y) « L{y\r]{e)) X 7t{0) (1) 
oc exp{-ia;2(y - rj{e)f} x eM-'^d^}- 

Thus an evaluation of the posterior requires a run of the forward model. While this simple, 1- 
dimensional density is trivial to evaluate, many inverse problems have to deal with a large model 
parameter vector (dimensions ranging from 10 to 10^) as well as a computationally demanding 
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forward model that may take a long time to evaluate (our experience ranges from seconds to 
weeks). 



1.1.1 Using a Gaussian process emulator 

While Markov chain Monte Carlo (MCMC) remains a popular approach for exploring the resulting 



posterior distribution ( Kaipio and Somersalo , 2004 ; Tarantola , 2005 ) , the demands required by the 



size of the model parameter vector and the computational demands of the forward model have 
inspired recent research focused on overcoming these hurdles. These research efforts range from 
response surface approximation of the forward model r/(-) ([Kennedy and O'Hagan 2001: Higdon 



et al. 



et al 



2005), to constructing reduced, or simplified forward models (Galbally et al 



2010), to polynomial chaos approximations of the prior model (Ghanem and Doostan 



et al., 2009). 



Marzouk and Najm 2009), to exploiting multiple model fidelities (Christen and Fox 



2010 



Lieberman 



2005 



2006 



Efendiev 



GP emulator-based solution 
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The forward model ry(-) is run at 4 values of the 
input parameter, producing model observations 
77° (black dots). A GP prior vr(77(-)) is specified 
for the forward model. 



likelihood (measurement) 
prior for r;(-) 
likelihood (simulations) 
prior for 6 
posterior 



Ad) 



Conditioning on both the physical observation 
(light, horizontal band) and the four model 
runs, the posterior distribution for both ?;(■) and 
9 is produced (shaded density). 



Figure 2: Using a Gaussian process prior for the forward model to reduce the number of model runs necessary for 
posterior exploration for the simple inverse problem. 



Gaussian processes (CPs) are commonly used to emulate the computer model response, pro- 
ducing a probabilistic description of the response at untried parameter settings - see [Kennedy and 
O'Hagan (2001) or Bayarri et al. (2007) for just a couple of examples. This basic approach is 
depicted in Figure [2[ In this case, a GP prior is used to model the unknown function 7]{-) and 
a collection of forward model runs r]° = (?/(6'°), . . . ,ri{6'^)y, over a collection of input parameter 
settings 6° = {91, . . . , 0^)' , are used to infer the forward model response at untried input parameter 
settings. 

Thus the basic fomulation ([T]) is augmented to incorporate this GP prior for the forward model 

GP(m(-),C(-,-)), 



where the mean function m(-) may be a constant (Sacks et aH 



or a more complicated regression function (Craig et al. 



2001 



1989 



Kennedy and O'Hagan 2001 ), 



Vernon et al. 2010), and the covari- 



ance function C(-,-) is typically of product form, requiring just a single additional parameter for 
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each dimension of the input parameter. For this simple problem, we take the mean and covariance 
functions as fixed, leading to the posterior form 



TT 



Here fig and vq are the mean and variance given by the GP model after conditioning on the forward 
model runs 7]° 

lie = C{e,e°)C{9°,9°)-\7]° -m{9°)) + m{e) 
VB = C{9,9) -C{9,9°)C{e°,9°)-^C{e°,9), 

where, for example, C{6°,9°) produces a m x m matrix whose ij entry is C{9°,9j) and C{9,9°) 
produces a m-vector whose jth element is C{9,9°). See Higdon et al. (2005) for details regarding 
the posterior specification when the GP (and other) parameters are not taken as fixed, and the 
input parameter is multivariate. 



1.1.2 Using the ensemble Kalman filter 

Below we briefly describe two basic variants of the the EnKF for computer model calibration, 
differing in how they use the ensemble of model runs to approximate, and represent, the result- 
ing posterior distribution. In both cases, an ensemble of draws 9° from the prior distribution of 
the model parameter are paired with the resulting simulation output to produce an ensemble of 
{9°,r](9°)) pairs, from which the sample covariance is used to produce an approximation to the 
posterior distribution. Hence we treat the input parameter settings 9° , ... ,9'^ as m draws from the 
prior distribution tt{9). Note that even though the distribution of the simulator response r]{9) is 
completely determined by the distribution for 9, the EnKF uses a joint normal model for (9,rj(9)) 
to motivate its calculations. 

Next we describe two variants of the EnKF algorithm for computer model calibration. One uses 
a Gaussian representation of the posterior distribution, the other uses an ensemble representation. 



Gaussian representation 

The first approach fits a multivariate normal distribution to the ensemble for (9° ,ri{9°)). The 
algorithm is depicted in the left frame of Figure [3] and described below. 

1. For each of the m simulations form the ensemble of joint vectors 

Jo)),A; = l,...,m. (2) 

With these m vectors, compute the sample mean vector /ipr and sample covariance matrix 
Spr. For the simple inverse problem here, //pr is a 2-vector and Spr is 2 x 2, but this recipe 
is quite general. 

2. In this simple inverse problem, the physical observation y corresponds to the 2nd element of 
the joint {9,r]{9)) vector. Take H be (0, 1)' to be the observation matrix. The likelihood can 
be written 

HvWe)) ^,J-\L-H (^^^ ) ' (, - H ) I (3) 
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More generally the observation operator H can select elements of r]{9) that are observed, or 
can be specified to interpolate between values of the simulator output. 

3. Combining the normal approximation to the prior with the normal likelihood results in an 
updated, or posterior, distribution for {9, rj) for which 

\y ^ ^(Mpost, Spost), (4) 

where 
and 

Aipost = Spost (SprVpr + H'T,~^y) . (6) 

Note that the posterior mean can be rewritten in a form more commonly used in Kalman 
filtering 

where T,p^H' {HT,p^H' + ^y)~^ is the Kalman gain matrix. 

The joint normal computations used here effectively assume a linear plus Gaussian noise relationship 
between 6 and r]{6), inducing a normal posterior for 6. 



EnKF Gaussian approximation EnKF ensembie approximation 




Figure 3: Left: Gaussian representation of the posterior distribution for {6,'q{6)) resulting from the ensemble 
Kalman filter (EnKF). The approximate normal prior distribution for {6,rj(6)) is depicted by the black ellipse, 
estimated from the ensemble (circle plotting symbols). The resulting posterior distribution is approximated as 
normal, depicted by the gray ellipse. The marginal posterior for 9 is given by the shaded density. Right: the 
ensemble representation of the posterior distribution for {9, r]{9)) resulting from the EnKF. Here the updated sample 
(gray dots) are approximate draws from the posterior distribution. 
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Ensemble representation 



The second approach is basically the usual EnKF as applied to time evolving systems, but here, 
only for a single time step. The goal is to perturb each member of the ensemble (^^j ^(^fc))) in order 
to produce an updated member {rj^\ G}^^) whose mean and variance match the posterior produced 
by the Gaussian representation EnKF described above. This updated member is not produced with 
the simulator so that rj^^ will not be equal to the simulator evaluated at updated parameter value 



Here we describe the perturbed data version of the EnKF described in 



Evensen 



V{0 

A number of variants of this basic approach exist; see Anderson (2001 ) and Szunyogh et a 
for example. The algorithm is given below. 

1. Construct the sample covariance matrix Spj- as in Step 1 of the previous algorithm. 

2. For k = 1, . . . ,m do: 

(a) Draw a perturbed data value ~ ^{Ui ^y)- 

(b) Produce the perturbed ensemble member 



(2009b) 



(2008), 



(1 



1) I = Spost ( ( ^^^^^ ) + H'^y'yk ) . 



(7) 



where Spr and Spost are defined in the previous algorithm. Note this perturbation of 
the ensemble member can be equivalently written using the more standard Kalman gain 
update: 

/ / ao \ 

(8) 




) + Ep,H'{H^p,H' + ^yr\yk - r,{9l)) 




3. Treat this updated, m member ensemble 

,{1 

1) I , /c = 1, . . . ,m. 

■ / 

as draws from the updated, posterior distribution for [9, rf) given the initial ensemble {9°,r]{9°)) 
and the physical observation y. 

This approach uses a Bayesian update of two normal forms, with each ensemble member updated 
separately. Here the normal prior is centered at the ensemble member, and the normal likelihood is 
centered at the perturbed data value, rather than at the ensemble mean and the actual data value. 
This update sets the new ensemble value to the mean of this resulting combination of 

normal distributions. 

This produces a posterior ensemble for the joint distribution of (9, t]), given by the gray dots in 
the right hand frame of Figure [3j Hence the difference between these two representations can be seen 
Figure [3] - compare the gray ellipse in the left frame, representing the updated normal posterior, 
to the gray dots in the right frame, representing draws using this ensemble representation. 

Note that if we take (0^, r/^) to be a draw from a distribution with mean fipj- and variance Spr, 

applying |7) - or equivalently (8) - produces a random variable {9^\'r]^^) with mean and variance 
Hence the mean and variance of the ensemble members {9^\'r]^^) matches 



given in mh and 



§, 
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that of the Gaussian representation of the EnKF Q. Even though the first and second moments of 
the two EnKF representations match in distribution, the ensemble representation appears to better 
capture the true, right skewed posterior (compare to Figure [T]). 

A two-stage approach 

Figure |4] shows how one can repeatedly apply the EnKF to improve the accuracy of the of the normal 
representation of r](9) where the posterior mass for 6 is concentrated. This iterative strategy is closer 
to the original use of the EnKF for state-space estimation in non-linear, dynamic systems. Also, 
this two stage approach easily generalizes to additional stages. 

For this two-stage EnKF, we artificially break the information from the likelihood into two even 
pieces 

L{y\riie)) oc exp - x exp - 7]{e)f 

as if y were observed twice, with twice the error variance. Then the EnKF is first applied to one of 
these y values, with twice the error varriance, producing an ensemble representation , . . • , Om^ of 
the posterior distribution for 9 given this partial piece of information. Next, the forward model is run 
again at each of these new parameter settings, producing the ensemble {9^\r][9^^)), k = 1, . . . ,m. 
This new ensemble is now the starting point for a second EnKF update, again using y with twice 
the error variance. 

This second update can produce a Gaussian representation (the gray ellipse in the right frame of 

I L (2) (2) 

Figure 4), or an ensemble representation [9j^ ,r]j^ ), k = 1, . . . ,m (the gray dots in the right frame 
of Figure |4]). As can be seen in the right frame of Figure |4| the second Gaussian representation of 
the relationship between 9 and r]{9) is more accurate because the {9^^\r](9^^^)) ensemble covers a 
narrower range, over which r]{9) is more nearly linear. 

Clearly, the choice of using two even splits of the likelihood information is somewhat arbitrary - 
both the number of splits and the partitioning of information to each split could be made in many 
ways. The cost of additional forward model evaluations has to be weighed against the benefits of 
a slightly more accurate Gaussian representation of r]{9) over a restricted range of values for 9. 

1.1.3 Embedding the EnKF into a Bayesian formulation 



As noted in a number of references (Anderson and Anderson 1999; Shumway and Stoffer, 2010 



Stroud et al. , 2010), the EnKF can be embedded in a likelihood or Bayesian formulation. For this 



simple inverse problem lends itself to the Bayesian formulation below, 

sampling model: y\rj{9) ~ N{r]{9),ay) 

prior model: (6*, 7?) ~ Af(/ipr, Spr), 

or, equivalent ly 

y\rj{9) ~ Nir]i9),al) 

r]\9 ~ N ^/ipr2 + ^pr22^pr2l(^ — A^prl), ^pt22 — 5]pr2lSpj.\ii;prl2 
9 ~ A^(/Xprl, Sprll). 
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1st stage EnKF ensemble approximation 2nd stage EnKF approximation 




-2-1012 -2-1012 

e e 



Figure 4: A two-stage EnKF solution to the simple inverse problem. Here the EnKF is applied twice, using the 
same observation y, but assuming it is observed with twice the variance. Left: at the first stage, the ensemble 
representation is used, assuming the observation has twice the variance, giving a data uncertainty that is a factor 
of \/2 larger than in the previous figures. Right: the second stage starts with the updated ensemble, evaluatingng 
the forward model at each 0^^', producing a new ensemble 77(6l<'^)), . . . , (sL^', r;(0ii>)). This new ensemble is 

updated once more, again using the data with twice the variance. The resulting uncertainty can be represented with 
a Gaussian distribution (shaded density) or an ensemble {9^^\ri^^), . . . , {9^\ri^) (shaded dots). 



In looking at the prior specification of 'q\9 above, it's apparent that the mean is just the Hnear 
regression estimate of r/ given 6. Hence where a GP model is used in the formulation described 
in Figure [2j the EnKF implicitly uses a linear regression-based emulator. While this simple form 
can only account for linear effects and no interactions, the tradeoff is that this emulator can be 
estimated quickly, can handle large ensemble sizes m, and can handle moderately high-dimensional 
input parameter, and output spaces. 

The EnKF uses the initial sample of model runs {9°,ri{9'^)),...,{6'^,r]{9'^)) to produce the 
standard plug-in estimates for /ipr and Spr - the sample mean and covariance. In static inverse 
problems, where quick turn-around of results isn't crucial, one could specify priors for these pa- 
rameters, producing a more fully Bayesian solution. An obvious choice might take vague, normal 
prior for /Xpr, and an inverse wishart for Spr (West and Harrison, 1997), if m is sufficiently large 
relative to the dimensionality of y and 9. 

In cases where the dimensionality of /ipr and Spr is large (much larger than the ensemble size 
m) covariance tapering, or some other form of localization is used to deal with spurious correlations 



produced in the standard sample covariance estimate ( Furrer and Bengtsson , 2007 ; Evensen 2009b 



Stroud et al. , 2010). The above specification suggests the use of variable selection (Wasserman 



2000 Tibshirani 1996), or compressed sensing (Baraniuk 2007) could make a viable alternative 



for estimating the regression function for rj given the ensemble draws for 9, producing the updated 
ensemble. Finally, we note that a bootstrap could be a useful tool for accounting for the uncertainty 
in the ensemble-based estimates for /Xpr and Spr since it does not require any additional model runs 
be carried out. 
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2 Applications 



This section describes three apphcations in the statistical analysis of computer models that make 
use of the EnKF. The first two are calibration examples, one taken from cosmology, the second from 
climate. The last explores how this EnKF representation can be used to for experimental design, 
determining optimal spatial locations at which to take ice sheet measurements. The goal of these 
examples are to suggest possible uses of EnKF ideas, rather than providing definitive analyses in 
problems involving inference with the aid of computationally demanding computer models. 



2.1 Calibration of cosmological parameters 

Perhaps the simplest cosmological model in agreement with available physical observations (e.g. 
the large scale structure of the universe, the cosmic microwave background) is the A-cold dark 
matter (ACDM) model. This model, controlled by a small number of parameters, determines the 
composition, expansion and fluctuations of the universe. 

This example focuses on model calibration, combining observations from the Sloan Digital Sky 



Survey ( Adelman-McCarthy et al. , 2006), giving a local spatial map of large galaxies, with large- 



scale A^-body simulations, controlled by five ACDM model parameters, evolving matter over a 
history that begins with the big bang, and ends at our current time, about 14 billion years later. 
An example of the physical observations produced by the SDSS are shown in the left frame of 
Figure [5] It shows a slice of the 3-d spatial map of large galaxies. Along with spatial position, the 
estimated mass for each of galaxy is also recorded. 

The computational model predicts the current spatial distribution of matter in the universe, 
given the parameters of the ACDM model, requiring substantial computing effort. For a given pa- 
rameter setting, a very large-scale A^-body simulation is carried out. The simulation initializes dark 
matter tracer particles according to the cosmic microwave background and then propagates them 
according to gravity and other forces up to the present time. The result of one such simulation is 
shown in the middle frame of Figure [sj Different cosmologies (i.e. cosmological parameter settings) 
yield simulations with different spatial structure. We would like to determine which cosmologies 
are consistent with physical observations of the SDSS given in the left frame of Figure [5j 




large scale structure data & simulations 




-3.0 -2.5 -2.0 -1.5 -1.0 -0.5 

wavenumber (loglO k) 



Figure 5: Left: Physical observations from the Sloan Digital Sky Survey (Credit: Sloan Digital Sky Survey). 
Middle: Simulation results from an A^-body simulation. Right: Power spectra for the Matter density fields. The gray 
lines are from 128 simulations; the black lines give spectrum estimates derived from the physical observations. 
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Table 1: ACDM parameters and their lower and upper bounds. 



parameter 



description 



lower upper 



n spectral index 0.8 1.4 

h Hubble constant 0.5 1.1 

as galaxy fluctuation amplitude 0.6 1.6 

^^CDM dark matter density 0.0 0.6 

r^B baryonic matter density 0.02 0.12 



Direct comparison between the simulation output and the SDSS data is not possible since the 
simulations evolve an idealized, periodic cube of particles corresponding to clusters of galaxies, 
while the SDSS data give a censored, local snapshot of the large scale structure of the universe. 
Instead, we summarize the simulation output and physical observations by their dark matter power 
spectra which describe the spatial distribution of matter density at a wide range of length scales. 

Computing the matter power spectrum is trivial for the simulation output since it is defined 
on a periodic, cubic lattice. In contrast, determining matter power spectrum from the SDSS data 
is a far more challenging task since one must account for the many difficulties that accompany 
observational data: nonstandard survey geometry, redshift space distortions, luminosity bias and 
noise, just to name a few. Because of these challenges, we use the published data and likelihood 



of [Tegmark et al. (2004) which is summarized by the black lines in the right hand frame of Figure 
[sj The resulting data correspond to 22 pairs {yi,ki) where yi is a binned estimate of the log 
of the power, and ki denotes the wavenumber corresponding to the estimate. The data vector 
y = (yi, . . . ,1/22)' has a diagonal covariance T,y. Two standard deviation error bars are shown in 
the right frame of Figure [5] for each observation. 



We take the ensemble produced in Heitmann et al. (2006) - a m = 128 run orthogonal array- 



based latin hypercube sample (LHS) over the 5-d rectangular parameter space detailed in Table 
[T] Since this sample was originally generated to produce a multivariate GP emulator - predicting 
the simulated matter power spectrum as a function of the 5-d parameter inputs - it is not a draw 
from a normal prior as is standard for EnKF applications. Nevertheless, this sample can be used 



to estimate //pr and Spr from Section 1.1.2 The restricted ranges of the parameters will need to 
be reconciled with the eventual normal description of the parameter posterior, or resulting EnKF 
sample. 

For each of the m = 128 parameter settings prescribed in the LHS, the simulation produces a 
55-vector of log power spectrum outputs, given by the gray lines in the right hand frame of Figure 
[5] Of the 55 elements in the simulation output vector, 22 of the elements are at the wavenumber 
k corresponding to the physical observations. Concatenating the parameter settings with the with 
the simulation output produces m = 128 vectors of length 5-1-55 

We take H to be the 22 x 60 incidence matrix, selecting the elements of the vector {9°,r]{9°)) 
that correspond to the physical observations. This, along with the physical observations y and 
corresponding measurement covariance T,y are the necessary inputs to carry out the Gaussian and 



ensemble representations of the EnKF described in Sections 1.1.2 
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Figure 6: The estimated posterior distribution for the model parameters (left) and the log power (right). Left: One- 
and two-dimensional marginals for the estimated posterior for the cosmological parameters. The upper triangle shows 
the 128 updated ensemble members (light circle plotting symbols) and a sample from the Gaussian representation 
(black dots). The lower triangle shows estimated 90% hpd contours for both the ensemble (light lines) and Gaussian 
representations (black lines). Right: posterior mean and pointwise 90% credible bands for log power spectrum for 
the matter density of the universe. Light lines give the estimate produced by the ensemble representation; black lines 
give the estimate produced by the Gaussian representation. 



The estimates of the posterior distribution for the 5-dimensional parameter vector is shown 
in the left frame of Figure [6| The presence of some shght skewness is noticible in the estimate 
produced by the ensemble representation. Also produced in these two estimation schemes is an 
estimate of the fitted log power spectrum, along with uncertainties, given in the right frame of 
Figure [6| These results, produced by the EnKF, can be compared to the posterior in Higdon et al. 
(2010), which was produced using a multivariate GP emulator, with a far more elaborate statistical 



formulation. The resulting posteriors are similar, but both EnKF estimates seem to "chop off" 
tails in the posterior for the cosmological parameters that are present in the GP emulator-based 
analysis. 



2.2 Optimal location of ice sheet measurements 

This second application comes from an ongoing effort to use the community ice sheet model (CISM) 



(Rutt et al. , 2009 Price et al. , 2011) along with physical measurements to better understand ice 



sheet behavior and its impact on climate. This study considers a model of an idealized ice sheet 
over a rectangular region which is flowing out to sea on one side, while accumulating ice from 
prescribed precipitation over a time of 1000 years. This implementation of the CISM depends on 
two parameters - 9i a constant in the Glen-Nye flow law (Greve and Blatter, 2009), controlling 



the deformation of the ice sheet, and 62 which controls the heat conductivity in the ice sheet. A 
few time snapshots of the model output are shown in Figure [7] for a particular choice of model 
parameters 61 and 6*2. 

While this configuration does not realistically represent important ice sheets in Greenland or 
Antarctica, it is a testbed where methodology can be evaluated for model calibration and/or planing 
measurement campaigns. After 1000 years, the thickness of the ice sheet could be measured to 
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200 years 500 years 1000 years 




Figure 7: Output from the idealized ice sheet model. The idealized ice sheet is described by height over a 36m x 
30m base, bounded on three sides by ledges. The fourth side is open to the ocean. Over the span of 1000 years, the 
ice flows into the ocean, while being replenished by a prescribed precipitation. Of interest in this application is the 
thickness of the ice sheet at 1000 years. 



inform about the model parameters 9i and ^2- The goal of this application is to use an ensemble 
of m = 20 model runs at different 6 = (^i, ^2) input settings to find a best set of 5 or 10 locations 
at which to measure the ice sheet thickness. 

The parameter settings and resulting ice sheet thickness (after 1000 years) for the m = 20 
model runs are shown in Figure [8j Thickness is produced on a 36 x 30 rectangular lattice of spatial 
locations. From this figure it's clear that the modeled ice sheet thickness is larger for smaller values 
of 9i , and larger values of 62 ■ 

The model runs produce an ensemble of 2 + 36 • 30 = p-vectors {OkjrjlOk)), k = 1, . . . ,m. We 
consider n ice thickness measurements taken at n of the 36 • 30 spatial grid locations given by the 
model. A given set of n measurement locations, is indexed by the n x p incidence matrix H, which 
will contain a single 1 in each of its n rows. Thus there are {p^2) possible measurement designs 
under consideration, each determined by which of the last p — 2 columns of H contain a 1. 

We use the Gaussian representation of the EnKF to describe the resulting uncertainty in 0, 
giving a simple means to compare designs which are determined by H. Assuming the n thickness 
measurements have independent measurement errors, with a standard deviation of one meter, means 
Tiy is the nxn identity matrix. Then the resulting posterior variance for the joint parameter-output 
vector is given by ([s]) 

^post = ^pr + ^J/ 

with the upper 2x2 submatrix of Spost describing the posterior variance for the parameter vector 
9. 

The sample covariance estimate for Spr, estimated from only m = 20 model runs, gives some 
spurious estimates for the elements of the covariance matrix, leading to aberrant behavior in esti- 
mates for conditional mean and variance for 9 given r]. If we define 

- (:) - ■ 

corresponding to the 2-vector 9 and the p — 2-vector r], a spatial tapering covariance matrix R{r) 
can be used to help stabilize these estimates (Kaufman et al. , 2008 Furrer and Bengtsson, 2007| > 
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Figure 8: Output from the ensemble of 20 model runs, showing thickness of the ice sheet after 1000 years, plotted 
in the (^i, S2)-parameter space. 6i controls deformation of the ice sheet; 62 controls heat conductivity in the ice sheet. 
Each image shows thickness as a function of spatial location; the center of the image marks the {9i, 62) input setting 
at which the model was run. The grayscale indicates log thickness. 



Hence we can produce an improved estimate 



T^r^'nir) = So R{r) 

where S denotes the sample covariance matrix from the samples rji, . . . ,rim, and o denotes the 
elementwise product of the matrix elements. Here we take R{r) to be the spatial correlation matrix 
induced by the isotropic exponential correlation function, with a correlation distance of r. The 
value for r is taken to be the maximizer of the likelihood of prior samples rji, . . . , r]m- 

1 f 1 1 

L(r) oc Yl \^vv(^)\~^ exp <^ -- {r]k - fJ-^)' - fJ-^) [ 

k=i ^ ^ 

Using this plug-in estimate for r, and treating Spr as known, the Bayesian D-optimal design 
that maximizes the prior-posterior gain Shannon information is simply the H that minimizes the 



determinant of the Spost ^ the upper 2x2 submatrix of Spost (Chaloner and Verdinelli, 1995). 
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Figure 9: Estimates of the Bayesian D-optimal designs for locations at which to measure the ice sheet depth for 
n = 3, 5, and 10. The estimates use a plug-in estimate for the spatial covariance distance of the covariance taper 
matrix, and Federov's exchange algorithm to carry out the optimization. 



Of course, since only a small number n of observations are likely to be taken, one need not 
compute using the full p x p matrix Spost- If we define Hrj to be the n x {p — 2) restriction of H, 
removing the first two columns of H, then the posterior covariance matrix for 9 can be written 



-"post 



^99 — '^erjH'{Hr^T.rj'q{r)H' + 



Here the computations require only the solve of a relatively small n x n system. 



We use the exchange algorithm of Fedorov (1972) to search for the design Hr^ that approxi- 
mately minimizes the determinant of Ep^g^ . The estimated optimal sampling locations for depth 
measurements with n = 3, 5, and 10, are shown in Figure[9j While the optimization algorithm only 
guarantees a local maximum, we tried a large number of restarts, with the configurations giving 
the minimal determinant of Sp^g^ shown in Figure joj 

2.3 Calibration of parameters in the Community Atmosphere Model 



This final application is an adaptation of the application described in Jackson et al. (2008), in which 



multiple very fast simulated annealing (MVFSA) was used to approximate the posterior distribution 
of climate model parameters. Here we use the EnKF to carry out model calibration, considering a 
more recent ensemble of 1,400 model runs, using the community atmosphere model CAM 3.1, as 
described in Jackson et al. (2008). In this application 15 model parameters are sampled uniformly 
over a 15-dimensional rectangle whose ranges are apparent in Figure 11 The model parameters, 



output fields, and corresponding physical observation fields are lis ted in Table [2} 

The computational model, described in detail in Jackson et al. (2008), produces a large number 



of outputs that could be compared to physical observations. We focus on a subset of the outputs 
(listed in Table [2]) explored in the original investigation. Each of these outputs is recorded as a field 
over the globe, averaged over 11 years (from 1990 to 2001), separately for each season (December 
- February, DJF; March - May, MAM; June - August, JJA; September - November, SON) . The 



images in Figure 10 show the two-meter air temperature observations. 

Rather that work directly with the model output and observed fields, we project these fields 
onto a precomputed empirical orthogonal function (EOF) basis, producing a small vector of weights 



one for each basis function - to represent each field (von Storch and Zwiers, 1999). As in Jackson 



et al. (2008), the EOF bases are computed from a long pilot run, separately from any of the model 



runs used to make the ensemble. The resulting weights for the two-meter air temperature are shown 



by the light/green dashes (model) and the black dots (observation) in Figure 10 We use 5 EOF 
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Table 2: Climate model inputs, outputs and physical data 



inputs 



outputs and physical data 



62 
O3 

04 

O5 
Oa 
Or 
0s 

eg 

ho 
hi 

9l2 
h3 
hi 

hs 



description 



effective radius of liquid cloud droplets over sea ice 

cloud particle number density over ocean & land 

effective radius of liquid cloud droplets over land 

time scale for consumption rate of deep CAPE 

cloud particle number density over warm land 

threshold for autoconversion of warm ice 

threshold for autoconversion of cold ice 

effective radius of liquid cloud droplets over ocean 

environmental air entrainment rate 

initial cloud downdraft mass flux 

low cloud relative humidity 

ice fall velocities 

low cloud relative humidity 

deep convection precipitation efficiency 

cloud particle number density over sea ice 



V4„y4 
V5,y5 
V6,y6 
V7,y7 



description 



shortwave cloud forcing 

precipitation over ocean 

two meter air temperature 

zonal winds at 300mb 

vertically averaged relative humidity 

air temperature 

latent heat flux over ocean 



basis elements for each output-season combination. Thus the model output ij and observation fields 
y are each summarized by a 7 x 4 x 5 vector of weights, corresponding to output, season, and EOF 
basis respectively. 

The long pilot run is also used to estimate the variation in the outputs expected just due to 
variation in climate. Thus for each output, season, and EOF, a variance c^iim ^^^^ estimated. 
We scale the EOF bases so that each cr^i^j^ is estimated to be 1. Thus, the error bars in Figure 10 



are ±2 because of this scaling. This scaling also makes the actual values of the y-axis in the figure 
essentially meaningless. 

Even with this variation estimated from the pilot run, it is expected that there will still be a 
discrepancy between the physical observations and the model output, even at the best parameter 
setting 9, for at least some of the outputs. Hence we specify to be the sum of the variance 
due to climate variation /140 and a diagonal covariance matrix that accounts for this additional 
discrepancy S^. For each output i we allow a different precision Aj for the discrepancy that is 
common across seasons and EOF bases. This gives 

T,s = diag (Aj^^ . . . , Af ^) ho 

so that 



The black dotted lines in Figure 10 show this additional uncertainty due to model discrepancy for 
the 2 meter air temperature, governed by A3. We specify independent T{a = 1,6= .001) priors for 
each Aj, z = 1, . . . , 7. 

In order to estimate these precision parameters, we note that the full 140-dimensional observa- 
tion vector y is modeled as the sum of normal terms 

y — 1] ~\~ Cclim ~l~ ^discrep 
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seasonal average 2 meter air temperature 1990 - 2001 




physical 
observation 
uncertainty due to 
weather variability 

uncertainty due to 
vi/eather variability 
and model 
discrepancy 



— prior ensemble 

— model output 



— posterior ensemble 

— model output 



Figure 10: Physical observations and uncertainty (black), prior simulations (light/green), and posterior predictions 
(dark/blue) for seasonal 2-meter air temperature, averaged from 1990 to 2001. The averages are computed for each 
season (DJF = Dec, Jan, Feb, and so on). The observed and simulated temperature fields are projected onto five 
EOF basis functions for each season, producing five EOF weights for each field. The basis was estimated using a 
single, long pilot run. The black dot shows EOF weights corresponding to the physical observations, the solid black 
line shows a 2-(Tciim bound for climate variation computed from the pilot run; the dashed black lines show additional 
uncertainty due to the estimated discrepancy error. The light/green lines are a sample of outputs from the ensemble 
of model runs. The dark/blue lines give the corresponding ensemble representation for the updated (i.e. posterior) 
model predictions. The scale of the y-axis has been standardized so that the estimated climate variance is one for 
each of the basis weights. The images above the plots show the physically observed two-meter air temperature fields. 



where rj ~ A^(/i^,S^^), with //^ and estimated from the prior ensemble as defined in 
eciim ~ N{0,Ii4o), and ediscrep ~ N{0,T,s). If we define 

V{X) = Y,r,r, + /l40 + ^5, 

we get the posterior distribution for the 7-vector A 

^(A|y) cc |y(A)riexp{-i(y-/z,)V(A)-i(y-/i??)} 

i=l 



We use the posterior mean as plug-in estimates for A, determining Sy. 

Now, given the dimension reduction from using the EOF bases estimated from the pilot run, 
the estimate for T,y, and the 1400 member ensemble of 15 + 140- vectors {9, rj{9)), and the ensemble- 
based estimates /ipr and Epr, the updated posterior distribution for r] and 9 is computed using 



the ensemble representation. The dark/blue dashes in Figure 10 show the posterior ensemble for 
the model outputs in the EOF weight space for the two meter air temperature. Figure 11 shows 
the posterior ensemble of parameter values 9. The prior ensemble was sampled uniformly over the 
15-dimensional rectangle depicted in the figure. 
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Figure 11: Ensemble representation for the posterior distribution of the 15 model parameters after conditioning on 
9 data fields. The parameter settings for the initial ensemble are uniform over the 15-dimensional rectangle depicted 
here. 



While the formulation presented here is very similar to that of Jackson et al. (2008), we used 
an additive discrepancy covariance matrix, with different precisions for each output type; theirs 
used a that is proportional to the estimated climate variation. Also, this analysis used fewer 
types of physical observations. The resulting posterior distribution for 9 is similar, with a bit more 
posterior spread in this analysis. Also, the EnKF analysis requires only an ensemble of model runs, 
with no need for the sequential sampling required for MVFSA. 

Finally, we point out that Annan et al. (2005) also use the EnKF to carry out parameter 
estimation on a climate model. That example uses a multi-stage estimation approach, collecting 
observations over ten successive years. That paper also uses synthetic observations so that Tiy can 
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be specified without the need for estimation. 



3 Discussion 



This paper highhghts a number of features of the EnKF from the perspective of model caUbration 
and shows examples of how it can be used in a variety of applications. Implicitly, the EnKF uses a 
multiple linear regression emulator to model the mapping between model parameters and outputs. 
This makes it easy for this approach to handle large ensembles ~ often a challenge for approaches 
that use GP-based emulators - as well as model outputs that are noisy or random. This also 
suggests regression-based approaches for dealing with high-dimensional input and output spaces 
may be helpful in EnKF applications. 

While the EnKF nominally starts with an initial ensemble from the prior distribution for 6, 
it's clear this prior will have little impact on the final results if the physical observations are fairly 
constraining, as in the examples presented here. The uniform designs used in the applications here 
have little impact on the posterior results. 

The results depend far more on specifications for covariance matrices, and how a large covariance 
matrix is estimated from a relatively small ensemble of model runs. We used likelihood and Bayesian 
approaches for estimation of covariance parameters; a variety of alternative approaches exist in the 



literature ( Tippett et al. 


2003 


Stroud and Bengtsson 


2007 


Evensen 


2009a 


Stroud et al. 


2010 


Kaufman et al. , 


2008 


)• 



The resulting posterior distribution for 6 tends to chop off tails that would be present in a more 
exact formulation. This is clear from comparing the analyses of the simple inverse problem laid 
out in Section |1.1[ This is largely due to the linearity of the regression based emulator implicitly 
used in the EnKF. We have also seen this phenomena in the ice sheet and cosmology applications 
when comparing to calibration analyses based on more exacting GP emulators. 

Finally we note that the ability of the EnKF to quickly provide "rough and ready" results makes 
it ideal for more computationally demanding tasks such experimental design or other optimization 
problems that require many iterations of the estimation process. The ice sheet application of 
Section [T2| is one such example. 
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